import numpy as np

PA = [0.3, 0.3, 0.2, 0.2]
PB = [0.4, 0.4, 0.1, 0.1]
PC = [0.2, 0.2, 0.3, 0.3]
PS = [[[0.2, 0.6, 0.2], [0.1, 0.3, 0.6], [0.05, 0.2, 0.75], [0.01, 0.1, 0.89]],
      [[0.6, 0.3, 0.1], [0.2, 0.6, 0.2], [0.1, 0.3, 0.6], [0.05, 0.2, 0.75]],
      [[0.75, 0.2, 0.05], [0.6, 0.3, 0.1], [0.2, 0.6, 0.2], [0.1, 0.3, 0.6]],
      [[0.89, 0.1, 0.01], [0.75, 0.2, 0.05], [0.6, 0.3, 0.1], [0.2, 0.6, 0.2]]]


def normal(res):
    # 计算res列表中所有元素的和
    sum_res = sum(res)
    # 将res列表中的每个元素除以sum_res，得到概率分布
    normalized_res = [res[i] / sum_res for i in range(len(res))]
    return normalized_res

#Gibbs采样
def Gibs():
    # n=20000
    n=4999
    res=[0,0,0]
    XA,XB,XC,sAB,sAC,sBC=0,0,1,0,1,1
    for k in range(n):
        _PA=normal([PA[i]*PS[i][XB][sAB]*PS[i][XC][sAC] for i in range(4)])
        XA=np.random.choice(4,p=_PA)
        _PB = normal([PB[i] * PS[XA][i][sAB] * PS[i][XC][sBC] for i in range(4)])
        XB = np.random.choice(4, p=_PB)
        _PC = normal([PC[i] * PS[XA][i][sAC] * PS[XB][i][sBC] for i in range(4)])
        XC = np.random.choice(4, p=_PC)
        sBC=np.random.choice(3,p=PS[XB][XC])
        res[sBC]+=1
    return normal(res)

#似然加权采样方法
def likehood_weighting():
    # n=20000
    n=5000
    res=[0,10,0]
    for i in range(n):
        w=1
        XA = np.random.choice(4, p=PA)
        XB = np.random.choice(4, p=PB)
        XC = np.random.choice(4, p=PC)
        w=w*PS[XA][XB][0]
        w=w*PS[XA][XC][1]
        sBC=np.random.choice(3,p=PS[XB][XC])
        res[sBC]+=w
    return normal(res)

#拒绝采样方法
def reject_sampling():
    # n=20000
    n=5000
    res=[0,0,0]
    for i in range(n):
        XA=np.random.choice(4,p=PA)
        XB=np.random.choice(4,p=PB)
        XC=np.random.choice(4,p=PC)
        sAB=np.random.choice(3,p=PS[XA][XB])
        sAC=np.random.choice(3,p=PS[XA][XC])
        sBC=np.random.choice(3,p=PS[XB][XC])
        if sAB==0 and sAC==1:
            res[sBC]+=1
    return normal(res)

#贝叶斯网络 精确算法,
def direct_cal():
    res=[0,0,0]
    for  XA in range(4):
        for XB in range(4):
            for XC in range(4):
                for sBC in range(3):
                    res[sBC] += PA[XA]*PB[XB]*PC[XC]*PS[XA][XB][0]*PS[XA][XC][1]*PS[XB][XC][sBC]
    return normal(res) #normal(X)=X/sum(X)将计数变成概率


if __name__ == '__main__':
    print("direct_cal algorithm#贝叶斯网络 精确算法:",direct_cal())
    print("reject_sampling algorithm#拒绝采样方法",reject_sampling())
    print("likehood_weighting algorithm#似然加权采样方法:",likehood_weighting())
    print("Gibs algorithm#Gibbs采样:",Gibs())